{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TOqzF3QC_p6U"
      },
      "source": [
        "# Tutorial: Analysis\n",
        "\n",
        "[![Open in Bohrium](https://cdn.dp.tech/bohrium/web/static/images/open-in-bohrium.svg)](https://bohrium.dp.tech/notebook/52a1b090ba2a4761836b5f19df1814e7)\n",
        "\n",
        "In this tutorial, we analyze a simple trajectroy as an example. The trajectory is taken from the real case [the thermal decomposition of CL-20](cl20.md).\n",
        "\n",
        "Download this trajectory using `wget`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8tAra7vQ_JOL",
        "outputId": "d44cc292-1e6c-4735-88e1-21c685b6d885"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "--2022-11-05 08:51:40--  https://raw.githubusercontent.com/tongzhugroup/TRAJREAX/f10a5c2cab77d3f3b659d9dd08256ae7b27c2820/cl20/cl20.lammpstrj\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.109.133, 185.199.110.133, 185.199.111.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.109.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 920127 (899K) [text/plain]\n",
            "Saving to: ‘cl20.lammpstrj’\n",
            "\n",
            "     0K .......... .......... .......... .......... ..........  5% 4.24M 0s\n",
            "    50K .......... .......... .......... .......... .......... 11% 5.15M 0s\n",
            "   100K .......... .......... .......... .......... .......... 16% 24.7M 0s\n",
            "   150K .......... .......... .......... .......... .......... 22% 20.2M 0s\n",
            "   200K .......... .......... .......... .......... .......... 27% 7.99M 0s\n",
            "   250K .......... .......... .......... .......... .......... 33% 42.7M 0s\n",
            "   300K .......... .......... .......... .......... .......... 38% 38.7M 0s\n",
            "   350K .......... .......... .......... .......... .......... 44% 40.9M 0s\n",
            "   400K .......... .......... .......... .......... .......... 50% 54.6M 0s\n",
            "   450K .......... .......... .......... .......... .......... 55% 58.8M 0s\n",
            "   500K .......... .......... .......... .......... .......... 61% 9.66M 0s\n",
            "   550K .......... .......... .......... .......... .......... 66% 45.7M 0s\n",
            "   600K .......... .......... .......... .......... .......... 72% 38.5M 0s\n",
            "   650K .......... .......... .......... .......... .......... 77% 80.2M 0s\n",
            "   700K .......... .......... .......... .......... .......... 83% 70.4M 0s\n",
            "   750K .......... .......... .......... .......... .......... 89% 68.7M 0s\n",
            "   800K .......... .......... .......... .......... .......... 94% 67.0M 0s\n",
            "   850K .......... .......... .......... .......... ........  100% 79.9M=0.05s\n",
            "\n",
            "2022-11-05 08:51:41 (18.5 MB/s) - ‘cl20.lammpstrj’ saved [920127/920127]\n",
            "\n"
          ]
        }
      ],
      "source": [
        "%%bash\n",
        "wget https://raw.githubusercontent.com/tongzhugroup/TRAJREAX/f10a5c2cab77d3f3b659d9dd08256ae7b27c2820/cl20/cl20.lammpstrj -O cl20.lammpstrj"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PmUV6VDyALqi"
      },
      "source": [
        "Then use `reacnetgenerator` to analyze it. Note the order of elements (`H C N O`) MUST map to those in the trajectory."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AgVDORot8sfl",
        "outputId": "bd357812-f3f7-4dee-e5cf-6a793b542ee6"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "2022-11-05 08:51:46,154 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: ReacNetGenerator: an automatic reaction network generator for reactive\n",
            "molecular dynamics simulation.\n",
            "\n",
            "Please cite: ReacNetGenerator: an automatic reaction network generator\n",
            "for reactive molecular dynamic simulations, Phys. Chem. Chem. Phys.,\n",
            "2020, 22 (2): 683-691, doi: 10.1039/C9CP05091D\n",
            "\n",
            "Jinzhe Zeng (jinzhe.zeng@ustc.edu.cn), Tong Zhu (tzhu@lps.ecnu.edu.cn)\n",
            "\n",
            "==================\n",
            "Features\n",
            "==================\n",
            "* Processing of MD trajectory containing atomic coordinates or bond orders\n",
            "* Hidden Markov Model (HMM) based noise filtering\n",
            "* Isomers identifying accoarding to SMILES\n",
            "* Generation of reaction network for visualization using force-directed\n",
            "  algorithm\n",
            "* Parallel computing\n",
            "\n",
            "==================\n",
            "Simple example\n",
            "==================\n",
            "ReacNetGenerator can process any kind of trajectory files containing \n",
            "atomic coordinates, e.g. a LAMMPS dump file prepared by running “dump 1\n",
            "all custom 100 dump.reaxc id type x y z” in LAMMPS:\n",
            "$ reacnetgenerator --type dump -i dump.reaxc -a C H O --nohmm\n",
            "where C, H, and O are atomic names in the input file. Analysis report\n",
            "will be generated automatically.\n",
            "\n",
            "Also, ReacNetGenerator can process files containing bond information, \n",
            "e.g. LAMMPS bond file:\n",
            "$ reacnetgenerator --type bond -i bonds.reaxc -a C H O --nohmm\n",
            "\n",
            "You can running the following script for help:\n",
            "$ reacnetgenerator -h\n",
            "\n",
            "2022-11-05 08:51:46,154 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Version: 1.6.6.dev24+g068dac4  Creation date: 2018-03-11\n",
            "\rRead bond information and Detect molecules: 0timestep [00:00, ?timestep/s]\rRead bond information and Detect molecules: 2timestep [00:01,  1.40timestep/s]\rRead bond information and Detect molecules: 101timestep [00:01, 70.84timestep/s]\n",
            "\rSave molecules:   0%|          | 0/11 [00:00<?, ?molecule/s]\rSave molecules: 100%|██████████| 11/11 [00:00<00:00, 3376.31molecule/s]\n",
            "2022-11-05 08:51:47,881 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 1: Done! Time consumed (s): 1.726 (Read bond information and detect molecules)\n",
            "2022-11-05 08:51:47,941 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 2: Done! Time consumed (s): 0.060 (Merge isomers)\n",
            "\rHMM filter: 0molecule [00:00, ?molecule/s]\rHMM filter: 11molecule [00:00, 2628.61molecule/s]\n",
            "2022-11-05 08:51:48,139 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 3: Done! Time consumed (s): 0.198 (HMM filter)\n",
            "\rIndentify isomers:   0%|          | 0/11 [00:00<?, ?molecule/s]\rIndentify isomers: 100%|██████████| 11/11 [00:00<00:00, 585.62molecule/s]\n",
            "\rAnalyze atoms:   0%|          | 0/11 [00:00<?, ?molecule/s]\rAnalyze atoms: 100%|██████████| 11/11 [00:00<00:00, 7775.08molecule/s]\n",
            "\rCollect reaction paths:   0%|          | 0/288 [00:00<?, ?atom/s]\rCollect reaction paths: 100%|██████████| 288/288 [00:00<00:00, 11114.42atom/s]\n",
            "\rAnalyze reactions (A+B->C+D):   0%|          | 0/100 [00:00<?, ?timestep/s]\rAnalyze reactions (A+B->C+D): 100%|██████████| 100/100 [00:00<00:00, 12790.24timestep/s]\n",
            "2022-11-05 08:51:48,583 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 4: Done! Time consumed (s): 0.444 (Indentify isomers and collect reaction paths)\n",
            "2022-11-05 08:51:48,651 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 5: Done! Time consumed (s): 0.068 (Reaction matrix generation)\n",
            "2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Species are:\n",
            "2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 1 [H]C12N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O\n",
            "2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 2 [H]C12NC3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O\n",
            "2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 3 O=NO\n",
            "2022-11-05 08:51:48,672 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: The position of the species in the network is:\n",
            "2022-11-05 08:51:48,672 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: {1: array([-5.95629110e-05, -1.00129074e-04]), 2: array([-0.59480018, -0.99989987]), 3: array([0.59485974, 1.        ])}\n",
            "2022-11-05 08:51:48,855 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 6: Done! Time consumed (s): 0.204 (Draw reaction network)\n",
            "2022-11-05 08:51:48,878 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Report is generated. Please see cl20.lammpstrj.html for more details.\n",
            "2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 7: Done! Time consumed (s): 0.087 (Generate analysis report)\n",
            "2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: ====== Summary ======\n",
            "2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Total time(s): 2.787 s\n"
          ]
        }
      ],
      "source": [
        "%%bash\n",
        "reacnetgenerator -i cl20.lammpstrj -a H C N O --type dump --nohmm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PO14aadxA4Cy"
      },
      "source": [
        "The results are shown in `cl20.lammpstrj.html`. In addition, we can view an overview of reactions in the trajectory in the command line:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GQNtq6mLAoJK",
        "outputId": "1fcc9399-5b5b-4fa6-b0aa-f753fcfb9fd2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1 [H]C12N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O->O=NO+[H]C12NC3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O\n"
          ]
        }
      ],
      "source": [
        "%%bash\n",
        "cat cl20.lammpstrj.reactionabcd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4xL87CdoBPQz"
      },
      "source": [
        "Here we can clearly see that the CL-20 molecule is decomposited into a nitrogen dioxide molecule.\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "authorship_tag": "ABX9TyN2IDYefbmLYArmNHl2104Q",
      "include_colab_link": true,
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.10.6 64-bit",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.10.6"
    },
    "vscode": {
      "interpreter": {
        "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
